Stable Skyrmions in two-component Bose-Einstein condensates 
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We show that stable Skyrmions exist in two-component atomic Bose-Einstein condensates, in the 
regime of phase separation. Using full three-dimensional simulations we find the stable Skyrmions 
with topological charges Q = 1 and 2, and compute their properties. With reference to these 
computations we suggest the salient features of an experimental setup in which they might realized. 



PACS Numbers: 03.75.Fi, 67.57.Fg 

Experimental advances in the formation and control 
of ultra-cold atomic gases are allowing detailed studies 
of the properties of Bose-Einstein condensates (BECs) 
containing multiple components. Multi-component con- 
densates have been formed by the simultaneous trapping 
and cooling of atoms in distinct hyperfine or spin levels 
and there are prospects of condensed mixtures of 
different bosonic atomic species |5|J| . The extra internal 
degrees of freedom introduced by the multiple compo- 
nents lead to a much richer phenomenology than in the 
single component case. 

One key feature of multi-component BECs is that, in 
general, their condensed groundstates are described by 
spinor order parameters ■ This raises the possibil- 
ity that these systems could support novel topological 
defects and textures. The existence of topological vor- 
tices |?],^| and monopoles JUS, in which the topology is 
imposed by boundary conditions, has been noted. How- 
ever, no localised topological textures, whose structure 
and topology is fixed only by simple energetic stability, 
have been identified. In Ref. jllj it was shown that a 
configuration with the topology of a "Skyrmion" [Q (a 
topological soliton of the S 3 — » S 3 map [|l3|) can be im- 
printed in a two-component Bose-Einstein condensate, 
using a carefully designed sequence of Rabi transitions. 
The subsequent time evolution indicated that this ob- 
ject is energetically unstable. In related work Q], it was 
demonstrated that a spherically-symmetric Skyrmion in 
a 2- or 3-component ferromagnetic BEC is also unstable 
to collapse. In Refs. a particular choice of inter- 

action coefficients was considered. However, the general 
multi-component system is described by all mutual two- 
body scattering lengths Given the range of param- 
eter space available it is natural to ask: Are there any 
conditions under which Skyrmions are stable? 

In this Letter, we present the results of extensive nu- 
merical studies of a fully 3-dimensional system that show 
that stable Skyrmions do exist in two-component BECs 
|L5j . We find that stable Skyrmions exist under the con- 
dition that phase separation occurs. We therefore expect 



Skyrmions to be stable in the experimentally relevant 
87 Rb \F = 2,m F = 1) |1,-1) system f|. As we will 
show, the Skyrmion can be viewed as a quantised vortex 
ring in one component close to whose core is confined the 
second component carrying quantised circulation around 
the ring. As such, the configurations closely resemble 
the "cosmic vortons" ]Tf^ | which may have formed in the 
early universe from superconducting cosmic strings fTij ] 
providing a further interesting link with between low- 
temperature laboratory experiments and cosmology (see, 
for example, Ref. Jl8[). 

We consider a two-component system, in which the 
number of atoms of each component is separately con- 
served on the timescale of interest, and assume that the 
many-body wavefunction can be written as a simple con- 
densate n 

N ± +N 2 

*({ri,r 2 ...})= J] [Vi(rO|l)+^(ri)|2)], (I) 

i=l 

where the states \a = 1,2) are the eigenstates of the two 
components (for simplicity of presentation, we take these 
to be two hyperfine levels). We are interested in stable 
stationary states, which are minima of the total energy 
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where the matrix U a p is determined by all mutual s-wave 
scattering amplitudes || . Since the number of atoms of 
each component is conserved, the minimisation is per- 
formed with constraints on 
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that is, having separate chemical potentials /ii ^ fi2- 
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The parameters entering (g) allow many different 
regimes. We will study solutions confined to a region 
close to the centre of a large trap that is loaded with 
component- 1 1). We may then neglect the effects of the 
confining potentials, and consider an infinite uniform sys- 
tem with boundary condition 
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where po is the density at the centre of the trap. Fur- 
thermore, we concentrate on situations in which the in- 
teractions are repulsive, and only weakly dependent on 
species Uu ~ U12 ~ 1/22- In this case, the stationary so- 
lutions vary on scales much larger than the healing length 
£ = (2mp Uh 2 )^ 1 ^ 2 . Consequently, we can neglect vari- 
ations in the total density p(r), and write 
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sin(0/2)e^ 



(5) 



in which we choose a parameterisation of the remaining 
freedom in terms of € [0, tt]; <j> a G [0, 2tt]. The field (§) 
and boundary condition (|J) allow a topological classifi- 
cation of (non-singular) field configurations, in terms of 
the winding number for the map S 3 — > S 3 
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We search for topological solitons by minimising, 
within each topological subspace Q, the energy functional 
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where A = -p 2 , [2f/ 12 — Uu — U22} with a constraint on 
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The constraint on N% is implicit in the choice of po- I n 
(R) the zero of the energy E is a constant which depends 
on Uu, U22 and the constrained numbers N\ and N2. 
There are two important length-scales in the problem 
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i?2 may be interpreted as the typical size of a cloud of 
component- 1 2); £a may be interpreted as the width of the 
transition region separating regions of component |1) and 
1 2) . Our neglect of the confining potential is justified pro- 
vided max(i?2j£A) is small compared to the overall size 
of the cloud; the restriction to constant density (|^) is jus- 
tified for £a 3> £j that is, |A| <C U. The dimensionless 
ratio T\ = (-R2 1 'Ca)j which is a measure of the amount 



of component- 1 2), is the parameter that characterises the 
solutions of (fjj), with the energies taking the form 



E Q (A,N 2 



h 2 p R 2 
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where £q(ji) is a dimensionless function. 

Consider first non-topological solutions (Q = 0). For 
the limit, V</>i = V02 = 0, the problem reduces to a 
3D double sine- Gordon model, which is known to have 
soliton solutions for A > 0, and r] > 1.92 (TJ). For 
large N2, the soliton can be viewed as a spherical do- 
main of component- 1 2) with radius ~ R2, separated from 
the region of component |1), by a domain wall of width 
~ £a- Therefore, it describes a spherical inclusion of 
component- 1 2) arising from phase- separation. 
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FIG. 1. The results of the 3D numerical relaxation for 
Q = 1 (on the left) and Q = 2 (on the right): (a) the 
iso-surfaces of the topological charge density (dark shading) 
and the number density of component- 1 2) (light shading); (b) 
the energy density on a plane containing the axis of symmetry; 
(c) the number density of component- [2) on the same plane 
and (d) the number density of component- 1 2), on a plane cut- 
ting the axis of symmetry. All figures correspond to n « 9. 

In contrast, there are no spherically symmetric topo- 
logical solutions (Q > 1); such solutions having at most 
axial symmetry. We have found such solutions for Q = 1 
and 2 using a full 3D simulation. Briefly, we evolved the 
full hyperbolic evolution equations of the system corre- 
sponding to the energy functional (Q) while constrain- 
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ing (||) using a Lagrange multiplier on a discretized grid 
with 100 3 points. Starting from a wide range of different 
initial conditions without any particular symmetry, but 
with the same value of Q, the kinetic energy was peri- 
odically removed so as to prevent oscillatory motion and 
effectively simulate gradient flow. This technique had 
previously been applied successfully to Skyrmions in the 
original model of Skyrme |2(J and also to the closely re- 
lated Hopf solitons |H] . The solutions evolved quickly to 
the required value of N2 and the subsequent relaxation 
was achieved after an acceptable period of time. As a 
simple test of the procedure we were able to reproduce 
quantitatively the results of Ref. jl9| for Q = 0. 

The results for Q = 1,2 are presented in Fig. [I] for 
the phase separation regime, A > 0, illustrating the ex- 
istence of stable topological solitons. We have plotted 
3D iso-surfaces of the topological charge density and the 
number density of species |2) (the integrands of (||) and 
(|J) respectively). Both illustrate the axial symmetry of 
the solutions: for this precise value of the threshold the 
topological charge density appears to be mostly localized 
around a central line in the shape of a "bolt" , although 
for a different threshold the figure would have appeared 
to be a vastly different since there is also some topological 
charge associated with the thin (~ £a) axially symmet- 
ric domain wall. The number density is clearly in the 
form of an axially symmetric ring encircling the "bolt" 
of topological charge density. 

Slices through the energy density (the integrand of (Q)) 
and the number density in the plane perpendicular to the 
ring illustrate the precise details of the solutions — in 
particular the two length-scales R2 and £a- The energy 
density, whose morphology is very similar to that of the 
topological charge density, is largely localized along the 
line through the centre of the vortex ring. However, one 
can clearly see that, emanating from this line, there is 
a shell of energy associated with the domain wall. The 
corresponding slice through the number density reaffirms 
this: it is almost zero away from the ring and rises sharply 
to close to one in the interior. An alternative view of the 
solution is a slice perpendicular to the line through the 
centre of the ring. 

In each case, the configuration can be viewed as a 
(unit) quantised vortex ring in component-|l), close to 
whose core is confined a circulating ring of component - 
|2) with quantised circulation (1 or 2 units, for Q = 1, 2). 
The kinetic energy of the vortex ring and the surface ten- 
sion of the domain wall separating the two components 
energetically favour shrinking the ring to zero size. This 
is balanced by the kinetic energy of the circulating core, 
allowing for stable configurations of non-zero radius. In 
this respect our solutions have similar qualitative form 
to the unstable configurations discussed in Ref. |11| for 
A = 0. Clearly, it is phase separation, induced when 
A > that allows stabilization of these solutions. 

Our 3D numerical studies indicate that the stable 



Skyrmions with Q = 1,2 are cylindrical^/ symmet- 
ric. To determine accurate values for the energies of 
these configurations, we have performed simulations us- 
ing a cylindrical ansatz consistent with the 3D results 
(and consistent with the variational equation for </>2): 
9(p, z), <f>i{p, z), 4>2 — m\ [(p, X, z) are cylindrical polar 
co-ordinates]. When m ^ 0, the remaining 2D problem 
for 9(p, z), </>i(p, z) has finite energy solutions provided 



= on p = and at p + z* 



this allows the con- 



figurations to be characterised by the winding number, 
n, of 4>i around the half-plane. Thus, there is a topologi- 
cal classification of the cylindrical ansatz in terms of the 
pair of integers (m,n), which describe the circulation of 
the vortex ring in component-|l) (n) and the circulation 
of component- 1 2) in the core (to). The winding number 
for the S 3 — ► S 3 map (||) is Q = nm. 

We have determined the soliton energies from numer- 
ical minimisations using a standard conjugant gradient 
routine. The results agree, to within numerical accuracy, 
with those of the 3D simulations. For m = we recover 
the non-topological soliton. We find stable topological 
solitons for all m > 0, but only for \n\ = 1. The ener- 
gies of the solitons £ m ,n=i(v) are plotted in Fig. 0. For 
Q > 3 (up to the largest we have studied) we find stable 
topological solitons within the cylindrical ansatz. These 
solutions may not be stable to non-axisymmetric pertur- 
bations which is currently under investigation. 
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FIG. 2. Scaled energies £ m n=1 of Skyrmions as a function 
of 77. The circulation of the vortex ring is n = 1; the circu- 
lation of the core is m = 1 ... 7 (from bottom to top), such 
that the topological charge is Q = m. The lines terminate 
below some critical value of n at which the solitons become 
unstable. Results are from calculations using the cylindrical 
ansatz on a 100 x 200 grid with R2 = 80 lattice spacings. 

As described above, the Skyrmions can be viewed 
as generalisations of quantised vortex rings (for exam- 
ple, in superfluid He) to a two-component condensate. 
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One distinction we wish to emphasise is that the stable 
Skyrmions of Figs, [j], || do not move in space, whereas 
conventional vortex rings move at constant (non-zero) 
velocity. We can find Skyrmion configurations that do 
propagate, by adding a constraint to fix the impulse 
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which is conserved by the time-dependent Gross- 
Pitaevskii equation. Under these (dissipationless) dy- 
namics, one can show that the configuration that min- 
imizes the energy at fixed P translates in space at a 
constant velocity v = dE/dP, where E is its energy. 
Dimensional analysis allows the energy of a state with 
constraints on Q, N 2 and P to be written as £(3(77, p) 
in Eq. (ji"o|) , where the new dimensionless parameter is 
p^\P\/(hpl /3 N 2 2 /3 ). 



Figure f (c,d) shows the distribution of component-|2) 
that would be observed experimentally for these config- 
urations. Clearly, the essential condition for the forma- 
tion of stable Skyrmions is that there is phase- separation 
between the two components. We therefore anticipate 
stable Skyrmions to exist for small deviations from the 
limit A <C U studied. Phase separation occurs provided 
H Ul 2 > U n U 22 (equivalent to A > for A < U). 
The observation of phase separation in the system 87 Rb 
\F = 2, mi? = f) |1, — f) 0, makes this a candidate sys- 
tem in which the stable Skyrmions may be realised. 

Our research is funded by EPSRC (PMS) and PPARC 
(RAB). The parallel computations were performed at 
the National Cosmology Supercomputing Centre in Cam- 
bridge. 




FIG. 3. Scaled energies of the stable Skyrmions with Q = 1 
(solid line, left axis) and Q = 2 (dashed line, right axis) as 
a function of scaled impulse, p, for n = 9. The minima cor- 
respond to the stationary Skyrmions. Calculations were per- 
formed using the cylindrical ansatz on a 100 x 200 grid with 
R 2 = 80. 

In Fig. |^ we show the dispersion relations of Skyrmions 
with Q = 1 and 2. The stationary configurations of 
Figs. |] and || correspond to the local minima (v = 
dE/dP — 0). Away from these minima, the velocities 
of the Skyrmions are non-zero. For p 1, the configu- 
rations resemble conventional vortex rings in which the 
radius of the ring is large compared to the core-size. 

In conclusion, we summarise the conditions for obser- 
vation of stable Skyrmions. In the limit that we have 
studied A <C U, we find solutions provided (i) A > 0, 
and (ii) N 2 is sufficiently large that 77 is above the termi- 
nation points of the curves in Fig. ^. Under these condi- 
tions, we expect that atomic Bose-Einstein condensates 
which have been imprinted with a topological configu- 
ration, in the manner of JnJ, will relax via dissipation 
to the stable Skyrmion configurations we describe here. 
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